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ADAPTIVITY AND BLOW-UP DETECTION FOR 
NONLINEAR EVOLUTION PROBLEMS 


ANDREA CANGIANI*, EMMANUIL H. GEORGOULIS*, IRENE KYZAt, AND STEPHEN 

METGALFE* 

Abstract. This work is concerned with the development of a space-time adaptive numerical 
method, based on a rigorous a posteriori error bound, for a semilinear convection-diffusion problem 
which may exhibit blow-up in finite time. More specifically, a posteriori error bounds are derived in 
the L^{H^)-type norm for a first order in time implicit-explicit (IMEX) interior penalty 

discontinuous Galerkin (dG) in space discretization of the problem, although the theory presented is 
directly applicable to the case of conforming finite element approximations in space. The choice of 
the discretization in time is made based on a careful analysis of adaptive time stepping methods for 
ODEs that exhibit finite time blow-up. The new adaptive algorithm is shown to accurately estimate 
the blow-up time of a number of problems, including one which exhibits regional blow-up. 

Key words, finite time blow-up; conditional a posteriori error estimates; IMEX method; dis¬ 
continuous Galerkin methods 


1. Introduction. The numerical approximation of blow-up phenomena in par- 
tial differential equations (PDEs) is a challenging problem due to the high spatial 
and temporal resolution needed close to the blow-up time. Numerical methods giv¬ 
ing good approximations close to the blow-up time include the rescaling algorithm of 
Berger and Kohn [SHU] and the MMPDE method |S||1S]- There is also work looking 
at the numerical approximation of blow-up in the nonlinear Schrodinger equation and 
its generalizations [IIIIIIIIISIIIIIISSI. Other numerical methods for approximating 
blow-up can be found for a variety of different nonlinear PDEs miaiisKisiiioiiiaiiQ] 
and ordinary differential equations (ODEs) |26l |29l [44] . Typically, these numerical 
methods rely on some form of theoretically justified rescaling but lack a rigorous 
justification as to whether the resulting numerical approximations are reasonable. In 
contrast, our approach is to perform numerical rescaling of a simple numerical scheme 
in an adaptive space-time setting driven by rigorous a posteriori error bounds. 

A posteriori error estimators for finite element discretizations of nonlinear para¬ 
bolic problems are available in the literature (e.g., IHEllTKIIllIlllMlIlIllielliTlIllI). 
However, the literature on a posteriori error control for parabolic equations that ex¬ 
hibit finite time blow-up is very limited; to the best of our knowledge, only in [33] 
do the authors provide rigorous a posteriori error bounds for such problems. Using a 
semigroup approach, the authors of |33] arrive to conditional a posteriori error esti¬ 
mates in the (L°° )-norm for first and second order temporal semi-discretizations 
of a semilinear parabolic equation with polynomial nonlinearity. Conditional a pos¬ 
teriori error estimates have been derived in earlier works for several types of PDEs, 
see, e.g., [HI11I3I1IMIIST]; the estimates are called conditional because they only 
hold under a computationally verifiable smallness condition. 

In this work, we derive a practical conditional a posteriori bound for a fully- 
discrete first order in time implicit-explicit (IMEX) interior penalty discontinuous 
Galerkin (dG) in space discretization of a non self-adjoint semilinear parabolic PDE 
with quadratic nonlinearity. The choice of an IMEX discretization and, in particular. 
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the explicit treatment of the nonlinearity, offers advantages in the context of finite 
time blow-up — this is highlighted below via the discretization of the related ODE 
problem with various time-stepping schemes. The choice of a dG method in space 
offers stability of the spatial operator in convection-dominated regimes on coarse 
meshes; we stress, however, that the theory presented below is directly applicable 
to the case of conforming finite element approximations in space. The conditional a 
posteriori error bounds are derived in the L°°(L^)-|-L^(iJ^)-type norm. The derivation 
is based on energy techniques combined with the Gagliardo-Nirenberg inequality while 
retaining the key idea introduced in [33] — judicious usage of Gronwall’s lemma. A 
key novelty of our approach is the use of a local-in-time continuation argument in 
conjunction with a space-time reconstruction. Global-in-time continuation arguments 
have been used to derive conditional a posteriori error estimates for finite element 
discretizations of PDEs with globally bounded solutions, cf. [S1IM1I3I1- A useful 
by-product of the local continuation argument used in this work is that it gives a 
natural stopping criterion for approach towards the blow-up time. The use of space- 
time reconstruction, introduced in |34[ 186] for conforming finite element methods and 
in [101128] for dG methods, allows for the derivation of a posteriori bounds in norms 
weaker than L?{H^) and offers great flexibility in treating general spatial operators 
and their respective discretizations. 

Furthermore, a space-time adaptive algorithm is proposed which uses the condi¬ 
tional a posteriori bound to control the time step lengths and the spatial mesh mod¬ 
ifications. The adaptive algorithm is a non-trivial modification of typical adaptive 
error control procedures for parabolic problems. In the proposed adaptive algorithm, 
the tolerances are adapted in the run up to blow-up time to allow for larger absolute 
error in an effort to balance the relative error of the approximation. The space-time 
adaptive algorithm is tested on three numerical experiments, two of which exhibit 
point blow-up and one which exhibits regional blow-up. Each time the algorithm 
appears to detect and converge to the blow-up time without surpassing it. 

The remainder of this work is structured as follows. In Section]^ we discuss the 
derivation of a posteriori bounds and algorithms for adaptivity for ODE problems 
whose solutions blow-up in finite time. Section [^ sets out the model problem and 
introduces some necessary notation while Section^ discusses the discretization of the 
problem. Within Section [^ the proof of the conditional a posteriori error bound is 
presented. An adaptive algorithm based on this a posteriori bound is proposed in 
Section [^ followed by a series of numerical experiments in Section]^ Finally, some 
conclusions are drawn in Section [8| 

2. Approximation of blow-up in ODEs. Before proceeding with the a poste¬ 
riori error analysis and adaptivity of the semilinear PDE, it is illuminating to consider 
the numerical approximation of blow-up in the context of ODEs. To this end, we first 
analyse the ODE initial value problem: find it : [0,T] —>■ K such that 



( 2 . 1 ) 


m(0) = Mo, 


with p > 2 a positive integer and coefficients ai > 0, i = 1,... ,p — 1 and Op > 0 so 


that the solution blows up in finite time. Let T* denote the blow-up time of (2.11 


and assume T < T*. For t < T*, u{t) is a differentiable function [27] . 
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Let 0 < < T, 0 < fc < iV be defined by := -L u- with t° := 0 and 

_ T for some time step lengths ti- > 0, k = 1,..., N with 
the following three different one step schemes to approximate ( |2.1[ ). We set '■= uq 
and, for k = 1,..., N, we seek such that 

- - - - =F{U'^-\U^), ( 2 . 2 ) 

'^k 

with F one of the following three classical approximations of /: 

Explicit Euler F{U'^-\U'^) = /(C/'^'"^), 

Implicit Euler F{U’^-\U'^) = f{U’^), (2.3) 

Improved Euler F(C/'=-\ U'^) = 1/2 (/(U'^-^) -f / {U^-^ + Tkf{U^~^))) . 

2.1. An a posteriori error estimate. We begin by defining {7 : [0, T] — M by 

C/(t) :=4-i(i)C/"-'+4(i)t^", (2.4) 


where denotes the standard linear Lagrange interpolation basis defined on 

the interval i.e., U is the continuous piecewise linear interpolant through 

the points (4,4^), k = 0,...,N. Hence, (2.21 can be equivalently written on each 
interval as 


dt ^ ^ 


(2.5) 


Therefore, on each interval (4 the error e := u — U satisfies the equation 


de 

dt 


= f{u) - F(C/'=-i, [/'=) = f{U) + f'(U)e F Y. - F{U’^-\ U^), (2.6) 


i=2 


with denoting the order j derivative of /. Thus, upon defining the residual 
r]k := fiU) — F{U^~^,U^), we obtain on each interval the primary error 

equation: 


de 

dt 


Vk 


f{U)t 


E 

f=2 


/a)([/) 


j! 


Gronwall’s inequality, therefore, implies that 


|e(t)| < Hk{t)Gk(j)k, 


where 


Hk{t) 



\k^Hu)\ 

j! 


U-1 



Gfc :=exp^^^ J/'(4)|ds^, 
(l)k ■■= |e(t''“^)| + [ |? 7 fc| ds. 


(2.7) 


( 2 . 8 ) 


(2.9) 







4 


A. CANGIANI, E. H. GEORGOULIS, I. KYZA AND S. METCALFE 


From (2.81, we derive an a posteriori bound by a local continuation argument. To this 
end, we define the set 

h ■= e [t'"~'^,t^] : ^max |e(s)| < 4Gfe(/>/c|, 


( 2 . 10 ) 


for some 4 > 1 to be chosen below; note that 4 C [t^~^,t^\ and that 4 is closed 
since the error function e is continuous. The main idea is to use the continuity of the 
error function e and to choose 4 implying that 4 = t*] for each k = 1,..., 4. 

Further, we will choose 6k to be a computable bound of Hk thereby arriving to an a 
posteriori bound. 

Theorem 2.1 (Conditional a posteriori estimate). Let u he the exact solution 
of (13, {4^}^o approximations produced by ( |2.2[ ) and U the piecewise linear 
interpolant (|2.4[). Then, for k = 1,..., 4, the following a posteriori estimate holds: 


max |e(t)| < SkGkfik, 


provided that 4 > 1 chosen so that 


Y^iSkGk^kY-^ 

i=2 


|/«(4(s))| 




j! 


ds - log(4) = 0. 


( 2 . 11 ) 


( 2 . 12 ) 


Proof. For A: = 1,..., 4, let 4 be as in ( |2.10[ ) where <5^ > 1 is chosen to satisfy 

|/«(C/(s))| 


exp ('^{SkGkffkY ^ [ 


J'- 


ds < (1 - a)4, 


(2.13) 


for some 0 < a < 1. The interval 4 is closed and non-empty since t^~^ G Ik', hence, 
it attains a maximum t^ := max4- Suppose that t^ < t^. In view of the definition 
of Hk, we have 

Hk{tt) < exp (^Y.(SkGkcl>ky-^ < (l-a)4 <4, (2.14) 


as 4 G 4- Application of (2.14) to (2.8) yields 

max |e(t)| < dkGk4>k- 


(2.15) 


This implies that t^ cannot be the maximal element of 4 “ a contradiction. Hence, 
t^ = t^ and, thus, 4 = \t^~^,t^\. Considering the case with equality in (2.13) and 
taking a —>■ 0, we arrive at ( |2.12 ) and the proof is complete. □ 

Choosing 4 > 1 satisfying (2.12) is equivalent to finding a root (ideally the 
smallest one) of 

i—2 


in the interval (1, -|-c»). This is only possible if the coefficients of a;4 j = 1,... ,p — 1, 
are “sufficiently small”, i.e., only provided that the time steps length Tk is small 
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enough. In this sense, (2.11) is a conditional a posteriori bound. In practice, condition 
(2.12) is implemented by applying Newton’s method. If for some k the time step 
length Tk is not small enough, Newton’s method does not converge and the procedure 
terminates (cf. Algorithms 1 and 2 below). With the aid of the next lemma, we state 
a precise condition on the time step lengths which indeed ensures that (2.12) has 
a root 5k > 1. 

Lemma 2.2. If 1 ~ log(2;) 'W'i-th Cj > 0, 

j = 1,... ,p, p € N has a root in (1, +oo). 

Proof. We begin by noting that s{x) is continuous in [l,-|-oo) and differentiable 
in (l,-|-oo) with s(l) = > 0 ^^.d limj,_>_|_oo s(a:) = -Loo. Therefore s{x) 

has a root in (1, -Loo) if and only if it attains a nonpositive minimum in this interval. 
Differentiating s{x) gives s'{x) = jCjX^~^—x~^. Since the coefficients Cj satisfy 

jCjC < 1, we observe that 


^< 0 . 

i=i 

Also, lima;_>_|_oo s'{x) > 0. Hence, there exists a critical point a;* G [e, +oo) satisfying 


X^jCjxi = 1. 


All that remains is to prove that s(a;*) < 0. Indeed, (2.16) leads to 


'^Cjxi < ^jCjxi = 1 < log(x*), 
i=i i=i 


(2.16) 


where the last inequality holds because a;* > e. From the above relation, we readily 
conclude that s(a:*) < 0 and the proof is complete. □ 

The above lemma gives a sufficient condition on when (2.12) can be satisfied. In 
particular, condition (2.12) can always be made to be satisfied provided that the time 
step length Tk is chosen such that 


P 

^1 {Gkhey-^ [ (17(s)) I ds < 1. 


Returning back to Theorem |2.1[ we note that this gives a recursive procedure for 
the estimation of the error on each subinterval Indeed, the term |e(t*’“^)| in 

(pk is estimated using the error estimator from the previous time step with e(0) = 0. 


2.2. Adaptivity. Based upon the a posteriori error estimator presented in The¬ 
orem |2.1[ we propose Algorithm 1 for advancing towards the blow-up time. 

Assuming that the adaptive algorithm outputs successfully at time T = T(tol, N) 
for a given tolerance tol and after total number of time steps N then we are interested 
in observing the order of convergence as T —>■ T* with respect to N. To this end, we 
define the function A(tol, A) := \T* — T(tol, A)| where T* is the blow-up time of 
(2.1) and we numerically investigate the rate r > 0 such that 


A(tol, A) oc A-F 


(2.17) 
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Algorithm 1 ODE Algorithm 1 
1; Input: /, F, uq, ti, tol. 

2; Compute from 17°. 

3: while / |77i| ds > tol do 

Jto 

4: Ti •<— Ti/2. 

5: Compute from 17°. 

6; end while 
7: Compute Si- 
8: Set k = 0. 

9: while S/c+i exists do 
10: k^k + 1. 

11: Tfc+1 = Tfc- 

12: Compute from [/'=. 

13: while / |77fc+i| ds > tol do 

Jt>‘ 

14: Tfe+i Tfe+i/2. 

15: Compute from 

16: end while 

17: Compute Sk+i- 

18: end while 
19: Output: k, 


One may initially expect that r would be equal to the order of the time-stepping 
scheme used. To gain insight into the rate convergence of A, we apply Algorithm 1 
to (2.11 with f{u) = for p = 2,3 and u(0) = 1 for each of the three time-stepping 
schemes ( 2.3| ). The computed rates of convergence r under Algorithm 1 are given in 
Table 12.1 ' 


Table 2.1: Algorithm 1 Results 


Method 

p = 2 

p = 3 

Implicit Euler 

r « 0.66 

r « 0.79 

Explicit Euler 

r « 1.35 

r « 1.60 

Improved Euler 

r « 1.2 

r « 1.48 


Somewhat surprisingly at first sight, the explicit Euler scheme performs signifi¬ 
cantly better than the implicit Euler scheme. This fact can be explained by looking 
back at the derivation of the error estimator. The explicit Euler scheme always un¬ 
derestimates the true solution u [H]. 

This, in turn, implies that is correcting for the fact that Gk+i is underesti¬ 
mating the true blow-up rate resulting in a tight a posteriori error bound and, thus, 
explaining the high rate of convergence of A. When using the implicit Euler method, 
on the other hand, Gk+i overestimates the true blow-up rate [44] thereby conferring 
no additional benefit. 

Note also that for both the implicit and improved Euler methods, the rate of 
convergence r is less than their formal orders of convergence, i.e., first and second 
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order, respectively. Moreover, one would expect a faster approach to the blow-up 
time using the second order improved Euler compared to the first order explicit Euler 
scheme. This unexpected behaviour is due to the way the tolerance is utilized in 
Algorithm 1. Indeed, Algorithm 1 aims to reduce the error under an absolute tolerance 
tol; this is the standard practice in adaptive algorithms applied to linear problems. 
In the context of blow-up problems, however, the presence of the growth factor Gk+i 
cannot be neglected; requiring the adaptivity to be driven by an absolute tolerance 
in the run up to the blow up time results in excessive over-refinement and, thus, loss 
of the expected rate of convergence. To address this issue, we propose Algorithm 2 
which increases tol proportionally to Gk+i allowing for control of the relative error 
(cf. line 19 in Algorithm 2). 


Algorithm 2 ODE Algorithm 2 
1; Input: /, F, uo, ti, tol. 

2: Compute from t/°. 

3; while / |? 7 i| ds > tol do 

Jto 

4: Ti •<— Ti/2. 

5: Compute from 1/°. 

6 : end while 
7; Compute Si- 
8: tol = Gi * tol. 

9: Set k = 0. 

10; while Sk+i exists do 
11; k ■(- k + 1. 

12; Tfc-i-i = Tfc. 

13; Compute C/'=+i from U^. 

.A+i 

14; while / |77fc_|_i| ds > tol do 

Jt>‘ 

15; Tfe+I ^ r/j+i/2. 

16; Compute 17^+^ from U^. 

17; end while 

18; Compute Sk+i- 

19; tol = Gk +1 * tol. 

20; end while 
21; Output: fc, 


Table 2.2: Algorithm 2 Results 


Method 

p = 2 

p = 3 

Implicit Euler 

r « 1.00 

r « 1.00 

Explicit Euler 

r « 1.45 

r « 1.43 

Improved Euler 

r « 2.03 

r « 2.03 


The rates of convergence r of A under Algorithm 2 are given in Table |2.2[ The 
theoretically conjectured orders of convergence for both the implicit and improved 
Euler schemes are recovered while the explicit Euler method still outperforms its 
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expected rate. In the case p = 3 (cubic nonlinearity) and for the explicit Euler 
method only, Algorithm 1 converges somewhat faster than Algorithm 2. The reason 
for this behaviour is unclear and requires further investigation. 

3. Model problem. Let C be the computational domain which is assumed 
to be a bounded polygon with Lipschitz boundary dil. We denote the standard Th¬ 
inner product on w C by and the standard L^-norm by Ij-Hcj; when uj = fl 

these will be abbreviated to (•, •) and Ij-lj, respectively. We shall also make use of the 
standard Sobolev spaces along with the standard notation U’{lo) = 

1 < p < oo; := A: > 0; and LAo(fl) denoting the subspace of H^{Vt) 

consisting of functions vanishing on the boundary d^l. For T > 0 and a real Banach 
space X with norm || • ||x, we define the spaces L^{0,T\X) consisting of all measurable 
functions v ■. [0, T] —>■ A for which 

for 1 < p < + 00 , 
for p = + 00 . 

We also define H^{Q,T,X) := {u C L^(0,T; A) : ut G L^(0,T; A)} and we denote by 
(7(0, T; A) the space of continuous functions z; : [0,T] —>• A such that 

lkllc(o.r;X) := max ||t>(t)||x < oo. 


IpIIlp(o,T;X) 


\ 1/p 




IMIlp( 0 ,T;A) := ess sup \\v(t)\\x < oo, 

0<t<T 


The model problem consists of finding m : x (0, T] —>■ K such that 


—- eAu + a • Vm + f(u) = 0 

at 

u = 0 

u(-,0) = Mo 


in O X (o,r], 

on dn X (0,T], 
in n. 


(3.1) 


for f[u) = /o — and where uq G iLo(^)j e > 0, a G [(7(0, T; W^’°°(n))]^ and 
/o G (7(0, T; L^(n)). For simplicity of the presentation only, we shall also assume that 
V • a = 0 although this is not an essential restriction to the analysis that follows. 

The weak form of ( |3.1| ) reads: find u G L^(0, T; i7g (fl)) (7 i7^(0, T; L^(n)) such 
that for almost every t G (0,T] we have 


du 


— ) + B{t;u,v) + {fit;u),v) =0 Vm G 


(3.2) 


where 


B{t-,u,v) := {eXu —au)-Xvdx. 

Jn 

Under the above assumptions and for any t G (0, T] the bilinear form B is coercive in 
and i7o(U), viz., B(t;v,v) > e||VM|p, for all v G 

4. Discretization. Consider a shape-regular mesh C = {K} of with K de¬ 
noting a generic element that is constructed via affine mappings Fk : K ^ K with 
non-singular Jacobian where K is the reference triangle or the reference square. The 
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mesh is allowed to contain a uniformly fixed number of regular hanging nodes per 
edge. On we define the finite element space 

V^(C) := {v € L^n) : v\k o Fk G rP(k), K G C}, (4.1) 

with V^{K) denoting the space of polynomials of total degree p or of degree p in 
each variable if K is the reference triangle or the reference square, respectively. In 
what follows, we shall often make use of the orthogonal L^-projection onto the finite 
element space V^, which we will denote by 11^. 

The set of all edges in the triangulation Q is denoted by 5(C) while 5*"*(C) C 5(C) 
stands for the subset of all interior edges. Given K £ ( and E G 5(C), we set 
Hk ■= diam(itr) and He ■= diam(£'), respectively; we also denote the outward unit 
normal to the boundary dK by n^- Given an edge E G 5“‘(C) shared by two elements 
K and K', a vector held v G [i5^/^(n)]^ and a scalar held v G i5^/^(0), we dehne 
jumps and averages of v and v across E by 

{v} := ^(''^IsnA + vIsnA'): M ■=v|_EnA ' + '^\EnK' ' ^K', 

•= 2 + I'lsnA')) H '■='^\EnK'^K + '>^\EnK''^K'■ 

If if C do., we set {v} := v, [v] := v • n, {?;} := v and [u] := vn with n denoting 
the outward unit normal to the boundary d^l. The inhow and outhow parts of the 
boundary dQ. at time t, respectively, are dehned by 

:= {x G (951 : a.{x,t) ■ n(a;) < 0}, (951*„( := {x G 951 : a.{x,t) ■ n(a;) > 0}. 


Similarly, the inhow and outhow parts of an element K at time t are dehned by 
(9i5*„ := {a; G dK : a{x,t) ■ riKix) < 0}, 955*^4 := {a; G dK : a{x,t) ■ nxix) > 0}. 


We consider an implicit-explicit (IMEX) space-time discretization of (3.21 con¬ 


sisting of implicit treatment for the linear convection-diffusion terms and explicit 
treatment for the nonlinear reaction term which was shown to be benehcial in Section 
For the spatial discretization, we use a standard (upwinded) interior penalty dis¬ 
continuous Galerkin method, detailed below, to ensure stability of the spatial operator 
in convection-dominated regimes. 

To this end, we consider a subdivision of [0, T] into time intervals of lengths 
Ti,... ,TAr such that T some integer N > 1 and we set := 0 and 

k := ^ = 1,... ,fV. Let C° denote an initial spatial mesh of 51 associated 

with the time = 0. To each time k, k = 1, ..., N, we associate the spatial mesh 
of 51 which is assumed to have been obtained from by local rehnement and/or 
coarsening. Each mesh is assigned the hnite element space := V/j(C*) given 


by (jg. For brevity, let := a(-, k) and ^ ■= /(’, Finally, for t G 

r(t) will denote the union of all edges in the coarsest common refinement u 
of and (k- 

The IMEX dG method then reads as follows. Let C/° be a projection of ug onto 
V°. For k = 1,..., n, find G such that 


ut - 


k-1 


Tk 


+ B{k-Ulvt) + = 0, (4.2) 
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for all S where 




KhiUlvl) 



ds 


- ^ {eWU^^}-[vl]+{eWvt}-[U^n]^s. 


We shall choose t/° as the orthogonal L^-projection of uq onto V°, that is 17° := n°uo, 
although other projections onto V° can also be used. In standard fashion, the penalty 
parameter, 7 , is chosen large enough so that the operator B + Kh is coercive on 
(see, e.g., [T7]). 

5. An a posteriori bound. In the context of the elliptic reconstruction frame¬ 
work [33113S] , we require an a posteriori error bound for a related stationary problem. 
To that end, we consider a generalization of the error bound introduced in |48] : the 
proof of such bound is completely analogous and is, therefore, omitted for brevity. 

Theorem 5.1. Given t C (0,r] and g G L^(n), let it® G Hq{Q) be the exact 
solution of the elliptic problem 

B{t]u%v) = (g,i;), 

for all V G 77(1(11) and let u'f G Yh such that 

B{t; ul, Vh) + Kh{ul, Vh) = {g, vh), 

for all Vh G Yh be its dG approximation. Then the following a posteriori error bound 
holds for any 0 v G 77,1 (fl); 


/ B{u^ - ul,v) 

V v^llVitll 


^ ^ ~ll5 + ~ ^ ^ ^^a||[V?1^]|||; 

NGC 



The symbols < and > used above and throughout the rest of this section denote 
inequalities true up to a constant independent of the data e, a, /, the exact and 
numerical solutions it, Uh, and the local mesh-sizes and time step lengths. 

Definition 5.2. We denote by G V* the unique solution of the problem 

B{t'^-Ulvi) + Kh{Ulvl) = yvtGYl 


For /c > 1, we observe that A^ = —11^/^“^ — (Uj^ — lP^Uli~^)/Tk from (4.2). 

Definition 5.3. We define the elliptic reconstruction G 77,1(11), fc = 0,..., A, 
to be the unique solution of the elliptic problem 


B{t'^;w\v) = {A\v) Vi;e 771(11). 


ut 


Crucially, the dG discretization of the elliptic problem in Def inition |5.3| is equal to 
and so B{t^\w^ — Ul^,v) can be estimated by Theorem 


5.1 


K 


At each time step fc, we decompose the dG solution Uf into a conforming part 
, G 77,1(11) nV^ and a non-conforming part U^^^G 


J\ such that = Ul 


jjk 
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Further, for t S ^we define Uh{t) to be the linear interpolant with respect to 
t of the values and viz., 

Uhit) -tk-imt" 

where, as before, {£k-i,^k} denotes the standard linear Lagrange interpolation basis 
defined on the interval We define Uh,c{t) and Uh,d{t) analogously. We then 

decompose the error e := u — Uh = ec — Uh,d with Cc '■= u — Uh,c and we denote the 
elliptic error hy 6^ •= — U^. 

Theorem 5.4. Given t G there exists a decomposition of Uh, as de¬ 

scribed above, such that the following bounds hold for each element K G IJ 

llVGh,dll],< h-^\\[Uh]\\l, 

EcKe 

\\Uh,d\\l< hE\\[UhWE, 

EcKe 

\\Uh,d\\L°°{K) ^ ll[^^]llL“(ifE)’ 

where Ke ■= {[J E : K Cl E Eg £{C,^ U denotes the edge patch of the 

element K — the union of all edges with a vertex on dK. 

Proof See |30j for the first two estimates and [16] for the final estimate. □ 

Lemma 5.5. Let t G then for any v G i4o(^) have 

+B{t;e,v)-\-{f{t-,u)-f{t-,Uh),v) = (^-f{t;Uh)- - B{t;Uh,v). 


Proof. This follows from (3.2). □ 


From Lemma 15.51 we obtain 
de 


,v] -G B{t-,e,v)if{t;u) - f(t;Uh),v) =-[ A +f 


dUh 


df- 

+ B{t^; v) - Bft- Uh, v) + B{e-Uf,v) + - f{t- Uh),v) 


dt 


(5.1) 


which, upon straightforward manipulation, gives 

Uh),v) = - + f'"~^ + 

+ ih-iB{t'^-^-,e'^-\v) + ikB{t^-,9\ v) - B{t- Uh, v) + 

+ Utv) + - fit- Uh) + 4-i(4l'= - v). 


In what follows, it will be convenient to define the a posteriori error estimator through 
three constituent terms rji, rjA, and rjB- The first part of the estimator is the initial 
condition estimator, rji, given by 


VI ■= 



+ 


E ’^eWp^We 

eg£(C°) 


1/2 


Both remaining parts, VA and vb, are the sum of a number of terms related to either 
a space or time discretisation error, identified by a subscript S or T, respectively. In 
this way, for t G VA is given by 


VA '■= ik-lVSi,k-l + (-kVSi.k + VS2,k + VTi,k, 
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where 


»7Si 


h\ 




fel||2 

h\\\E 


agc'" 


»?S2.fe := 


while ? 7 _b is given by 


where 


E 


[U^k]\ 

11 + 

^||[a'= 44 || 

Bg£(C* 





E 

t^K 

f 

11 /'=- 

1 ^ 

^ agc'““^uc 

k 



= £- 1 /^ 

114-1 

(a'=- 

i-a 

)t^r '+40 


£;g£‘"‘(C'‘) 

1/2 


1/2 


iB ■= VSi.k + r]S4,k + 'nT 2 ,k, 


'ns3,k ■= 


E E ^^k^eWP] 

AgC'‘-iuC'' E( 1 Ke 


h\\\% 


1/2 


1/2 


with 


E hE\mt-Ut^)/ru]\\l 

^ Acr(t) 

:= 11/"-' - fit -. Uh ) + 4-1(21'= - ^'=-1) 

(Tif := 2||C/?,||ioo(^) + ||[17;i]||ioo(//g). 


With the above notation at hand, we go back to (5.2) and bound the first term 
on the right-hand side using the definition of , the orthogonality property of the 
L^-projection and the Cauchy-Schwarz inequality: 




dUh 

dt 


;< 77S2,fcv4||Vu||. 




Tk 


,u-n'=u 


(5.3) 


The next two terms give rise to parts of the space estimator via Theorem |5.1| 

tk-iB{t^-^-e^-\v)+ikB{t^-e\v) < (4-ir/Si.fc-i +4^751, fe)V^||Vu||. (5.4) 

Using the definition of the bilinear form B and the Cauchy-Schwarz inequality, the 
final four terms give rise to the time estimator: 


4-iS(4-i;[/ri,u)+4B(t'';C/^u)-B(t;C/^,u) <7 ?TiAv4||Vu||, 

(/'=-!-/(t; 44 + 4 -i(dl'= - A'=-i), 4 < r?T 2 ,fe|kl|. 


(5.5) 


Setting u = Cc in (5.21, using the results above along with the coercivity of the bilinear 
form B and the Cauchy-Schwarz inequality, we obtain 


1 d 

2 dt 


edl^ -b ellVedP -b - f{t;Uh),ec) < 


.dU, 


h,d I 


dt 


^T2.fc 


(4-i?7Si,fc-i + 4??Si,fc + 'riS2,k + ?7ri.fc)v4||Vec|| -b B{t; Uh,d, Sc)- 


(5.6) 
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Application of Theorem |5.4| implies that 

+ £||Vecf + {f{t;u) - f{t;Uh),ec) < + VT 2 ,k)\\ec\\ 

+ iik-ir]Si,k-l + (-kVSi,k + r]S2,k + i?Ti,fc)-\/e|| VCcH . 


Thus, we conclude that 


1 d 

2 dt 




1 

2 


r]A + VB\\ec\\- 


(5.8) 


We must now deal with the nonlinear term on the left-hand side of (5.8). We begin 
by noting that 


u) - fit] Uh),ec) = {f{t; Cc - Uh,d + Uh) - f{t\ Uh), e^) = Ti + T2, 


(5.9) 


where 

Ti := {2UhUh,d,ec) - iU^^^,ec), 

T 2 ■= —{2Uhec, Cc) + {2ecUh,d, Sc) ~ (cci Cc)- 
Upon writing the contributions to Ti elementwise and using Theorem |5.4[ we have 

/ \ 1/2 

|Ti|<( ^ {nUhWL^m + Wu^AL^mYwUhAK) l|ec|| 


(5.10) 


< 


^Sa.fcllecll- 

To bound T 2 , we use Holder’s inequality along with Theorem |5.4| to conclude that 

lAl ^ (2||U/j||ioo(n) + l![t^/i]||L=o(r(t)))||ec|P + ||ec||i3(n)- (5-11) 

Combining (5.8), ( |5.9| ), ( |5.10[ ) and (5.11) we obtain 


-||e,f + £||Ve,f < CAa + 2Cr^B\A\\ + 2an||e,f + 2||e,||i3(a), 


(5.12) 


with tro := 2\\Uh\\L°°{a) +C'l|[C^/i]||L“(r(t)) where C > 0 is a constant that is indepen¬ 
dent of e, a, /, u, Uh and the local mesh-sizes and time step lengths. For v G HA^), 
the Gagliardo-Nirenberg inequality < C'gn||i'||^|| V?;|| implies that 

||ee||i3(0) < CcNlle.rllVe.ll < |||Ve,f + (5.13) 

Thus, 

^||ec|P < -I-2C77_B||ec|| + 2crQ||ec|p-I-C gnS ^l|ec||^- (5-14) 

To deal with the L^-norms of Cc appearing on the right-hand side, we use a variant 
of Gronwall’s inequality. 

Theorem 5.6. Let T > 0 and suppose that Cq is a constant, Ci,C 2 € L^(0,T) are 
non-negative functions and that u G lU^’^(0,r) is a non-negative function satisfying 


A{T)<cI-\- f ci(s)m(s) ds-I- f C2{s)A{s)ds, 

Jo Jo 
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then 


u{T) < (^\co\ + ^ ci(s)ds^ exp Qy C 2 (s)ds 
Proof. See Theorem 21 i n |19l . □ 

Application of Theorem 5.6 to (5.14) for t G (t^~^,t^j yields 


|ec(t)|| < 'Hk{t)Gk^k, 


(5.15) 


where 


:= 


|ec(t 


k-l\ 


C 




Ta ds 


1/2 


C 


1 


r]B ds, 


Gk := exp 


TO ds , 


'Wfe(t) := exp ^ lledpds 


To remove the non-computable term Hk from (5.15), we use a continuation argument. 
We define the set 

Ik ■= {t e : ||ec||Loo((fe-i < SkGk^k}, 

where, analagous to the ODE case, 6k > I should be chosen as small as possible. Ik 
is non-empty (since G Ik) and bounded and, thus, attains some maximum value. 
Let t* = maxXfc and assume that t* < t*. Then, from (5.15), we have 

l|ec||L“(t''-i.f;L2(n)) < 'H{t*)Gk^k 

< exp (C^N^”^T-fe||edlioe(t^-i_t,.i2(o)))Sfe4>fc (5.16) 

< exp {C^Ne-'^Tk6lGl<i>l)Gk<^k- 
Now, suppose <5fc > 1 is chosen such that 

exp {C^N^-^Tk6lGl^l) < (1 - a)Sk, (5.17) 

for some 0 < a < 1 then ( |5.16[ ) gives 

(1 — a)SkGk^k < 6kGk^k, (5.18) 

which, in turn, implies that t* cannot be the maximal value of t in — a contradiction. 
Hence Ik = and we have the desired error bound once Sk is selected. Taking 

a —>■ 0, we can select J/c > 1 to be the smallest root of 

CGN£~^Tk6lGl<^>l - log(4) = 0. (5.19) 

Finally, we estimate $i. Application of Theorem |5.4| and the triangle inequality yields 

\\eMf<\mr + \\UH,Mr<cvi ( 5 . 20 ) 

Therefore, if we redefine $i to be 
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we have 


||ec(t^)|| < ||ec||L“(tO,ti;_L2(Q)) < 'Pi, 

where 'I'l := SiQi^i. In the same way, if we redefine 

1/2 


(5.21) 


■■= ('i’l-i+C [ f Vsds, 


d’k '■= SkGk^k, 


we have 


|ec(i*^)|| < ||ec||L“(i''-i,t'=;L2(o)) < 'pfc. 


(5.22) 


Hence, we have shown the following result. 


Theorem 5.7. The error of the IMEX dG discretization of problem (3.2), given 

1/2 


by (4.2), satisfies 


|e||L“(0,T;L2(Q)) < Tesssup f V /iB II [Ha] III) , 

o<t<T \ / 


Bcr(t) 


providing that the solution to (5.19) exists for all time steps. 

Proof. Follows from (5.22), the triangle inequality, and Theorem 5.4 □ 

The estimator produced above is suboptimal with respect to the mesh-size as it 
is only spatially optimal in the L^(77^)-norm. It is possible to conduct a continuation 
argument for the L^(i7^)-norm rather than the L°°(L^)-norm if one desires a spa¬ 
tially optimal error estimator; this is stated for completeness in the theorem below. 
However, the resulting S equation was observed to be more restrictive with regards 
to how quickly the blow-up time is approached. For this reason, we opt to use the a 
posteriori error estimator of Theorem 5.7 in the adaptive algorithm introduced in the 
next section. 


Theorem 5.8. The error of the IMEX dG discretization of problem (3.2), given 


by (4.2), satisfies 


rT X 1/2 Y . 

e(r)|p-f/ e||Vec|pdt) < ^ 4'^-f esssup I ^ /ie||[Ha]||| 

Jo J tx 0<t<T V 


1/2 


Furthermore, close to the blow-up time where ||e(r)|| = ||e||Bcx,(o_T;L 2 (n)) we have 

/ rT X 1/2 Y . X 1/2 

l|e||i~(o,T;L2(n)) + / e\\Vec\\^dt] <^Tfc-fesssup ^ /ib||[Ha]||| , 

where 'Pfc, fc = 1,..., TV, is defined recursively with dtg = Crji and 

1/2 


^k +C rjj^ds C 


/ \ 

Gk := exp(Tfe/2) exp ^ an dsj , 

4'fc := SkGk^k, 


11B ds 
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provided that <5^ > 1 which is the smallest root of the equation 

- log( 4 ) = 0 , 

exists for all time steps. 

Proof. The proof is completely analogous to that of Theorem |5.7| and follows from 
(5.12) by conducting a continuation argument for the )-norm. □ 


6. An adaptive algorithm. We shall now proceed by stating our space-time 
adaptive algorithm for problems with finite time blow-up. 

The algorithm is based on Algorithm 2 from Section 2 and the space-time adaptive 
algorithm for linear evolution problems presented in m- It makes use of different 


terms in the a posteriori bound in Theorem 5.7 to take automatic decisions on space¬ 


time refinement and coarsening. The pseudocode describing the adaptive algorithm 
is given in Algorithm 3. 

The term rjSi.k drives both local mesh refinement and coarsening. The elements 
are refined, coarsened or left unchanged depending on two spatial thresholds stol+ 
and stol~. Similarly, the term r]T 2 ,k is used to drive temporal refinement and coars¬ 
ening subject to two temporal thresholds ttol+ and ttol“. 

7. Numerical Experiments. We shall investigate numerically the a posteri¬ 


ori bound presented in Theorem 5.7 and the performance of the adaptive algorithm 
through an implementation based on the deal. II finite element library [4]. All the 
numerical experiments have been performed using the high performance computing 
facility ALICE at the University of Leicester. The following settings are common to 
all the numerical experiments presented. We use polynomials of degree five, hence 
p = 5 throughout. This particular choice provides a good compromise between run 
time and spatial discretization error for the problems considered. Furthermore, it 
permits us to analyse the asymptotic temporal behaviour of the adaptive algorithm. 

Correspondingly, we set 7 = 30 to ensure coercivity of the discrete bilinear form. 
We also set ttol“ = 0.01 * ttol+ and stol“ = 10“® * stol+ as, respectively, the 
temporal and spatial coarsening parameters. The initial mesh is chosen to be 
a 4 X 4 uniform quadrilateral mesh and the initial time step length ti is chosen 
so that the first computed numerical approximation is before the expected blow-up 
time. The unknown constants in the a posteriori bound are set equal to one as is the 


constant Cqn in (5.19); the above conventions are deemed reasonable for the practical 


implementation of the a posteriori bound. 

7.1. Example 1. We begin by considering a standard reaction-diffusion semi- 
linear PDE problem whose blow-up behaviour is theoretically well understood. This 
is given by setting U = (—4,4)^, e = 1, a = (0,0)^, /o = 0 and uq = 

The initial condition uq is chosen to be a Gaussian blob centred on the origin that is 
chosen large enough so that the solution exhibits blow-up; the blow-up set consists of 
a single point corresponding to the centre of the Gaussian. 

To assess the asymptotic behaviour of the error estimator, we fix a very small 
spatial threshold so as to render the spatial contribution to both the error and the 
estimator negligible. We then vary the temporal threshold and record how far the 
algorithm is able to advance towards the blow-up time. The results are given in Table 

rm 

For the present case (problems without convection), it is known that the solution 


to (3.2) has the same asymptotic behaviour as the solution to (2.1) with respect to the 


time variable m- Thus, we would expect an effective estimator to yield similar rates 
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Algorithm 3 Space-time adaptivity 
1 ; Input: e, a, /o, uq, 7 , ttol+, ttol“, stol+, stol^. 

2 : Compute C/°. 

3: Compute from Uj^. 

4: while / rj^ ^ds> ttol"*" OR max? 7 g ■^\k > stol+ do 
JtO K 

5: Modify C° by refining all elements such that rjg^ ilif > stol+ and coarsening 

all elements such that rjg i\k < stol“. 
ft" 

6 : if / rj^ ids > ttol^ then 

JtO 

7: Tl <r- ri/2. 

8 : end if 

9: Compute U^. 

10 : Compute Ul from U°. 

11 ; end while 
12 : Compute 

13: Multiply ttol+, ttol^, stol+, stol“ by the factor Gi. 

14: Set j = 0, C = C°- 
15: while 5j+i exists do 
16: j + 1 . 

17: T,+i=T,. 

18: Compute Ui'^^ from 

/.U+i 

19: if / j_i_i ds > ttol'*' then 

JtJ 

20: Tj + i ^ Tj + i/2. 

21 : Compute from 

22 : end if 

23: if / 77^2 j_i_i ds < ttol then 

Jv 

24: Tj+i ^ 2 Tj+i. 

25: Compute from 

26: end if 

27; Form from (ti by refining all elements such that j+iIk > stol+ and 

coarsening all elements such that < stol“. 

28: Compute from Ul- 

29: Compute 6j+i. 

30: Multiply ttol+, ttol^, stol+, stol“ by the factor Gj+i. 

31; end while 

32: Output: j,tG \\Uh{P)\\L’^(^ny 


for A, the difference between the true and numerical blow-up time, to those seen in 
Section Although the true blow-up time for this problem is not known, we observe 
from Table iTdl that 


l!t^7i||L“(0,r;L“(O)) OC 


From m, we know the relationship between the magnitude of the exact solution in the 
L°°(L°°)-noTT[i and the distance from the blow-up time. Thus, under the assumption 
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Table 7.1: Example 1 Results 


ttol+ 

Time Steps 

Estimator 

Final Time 

\\Uh{T)\\L^^n) 

1 

3 

9.5 

0.09375 

12.244 

0.125 

8 

24.6 

0.12500 

14.742 

(0.125)2 

19 

54.0 

0.14844 

18.556 

(0.125)3 

42 

66.7 

0.16406 

23.468 

(0.125)4 

92 

218.5 

0.17969 

32.108 

(0.125)3 

195 

1142.4 

0.19043 

44.217 

(0.125)3 

405 

1506.0 

0.19775 

60.493 

(0.125f 

832 

1754.1 

0.20313 

83.315 

(0.125)3 

1698 

5554.2 

0.20728 

117.780 

(0.125)9 

3443 

6020.4 

0.21014 

165.833 

(0.125)4° 

6956 

33426.7 

0.21228 

238.705 

(0.125)44 

14008 

36375.0 

0.21375 

343.078 

(0.125)42 

28151 

66012.8 

0.21478 

496.885 

(0.125)43 

56489 

157300.0 

0.21549 

722.884 


that the numerical solution is scaling like the exact solution, we have 
A(ttol+, A^) « ||w||^oo(g_y.^cx=(Q)) ~ l|t^/i|lL=°(o,T;L“(n))' 

Therefore, we conjecture that 

A(ttol"'", IV) cx 

Note that the conjectured convergence rate is slower than the comparable results in 
Section a possible explanation for this will be given in the concluding remarks. 

Next we investigate the numerical blow-up rate of ||t7/t(t)||L~(n). In particular, we 
are interested in checking if the numerical blow-up rate coincides with the theoretical 
one. For this particular example, it is well known, cf. [351 ISS]) that close to the 
blow-up time ||u(t)||icxj(f 2 ) behaves as 

ll'«(t)||L~(n) ^ T* -t ' 

where T* denotes the blow-up time. Let us denote the numerical blow-up time by 
t* which we compute as follows. For the last numerical experiment (with ttol'^ = 
(0.125)^^), we assume that there exists a constant Cn such that 

\\UHmL^m = CN^^, t = t^-\T. 

Then t* is computed by 

. T\\Uh{T)\\L^ia) - t^-^\\Uh{t^-^)\\L-(n) 

\\UhiT)\\La<>(^Q) - ||C/?,(t^“l)||Loo(Q) 

For this example, the above relation gives t* = 0.217055. 

To compute the numerical blow-up rate, we consider all the time nodes t^, k = 
0,... ,N with t^ = 0.21375 (corresponding to the numerical experiment with ttol“'" = 
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(0.125)^^). Then for every two consecutive times we assume that there exists 

a constant Ck such that 





1 

- t)P>‘ ’ 




and hence we compute pk by 

_ log(ll^hft'")IU°°(a)/l|b^/i(4-i)||L°°(n)) 

— t'^)) 


This produces a sequence {pk}k=i numerical blow-up rates. Since for this example 
the theoretical blow-up rate is one, for a “correct” asymptotic blow-up rate of the 
numerical approximation we expect pk to tend to a number close to one as k ^ N. 


This is indeed the case, as observed in Figure 7.1 


7.2. Example 2. Let fl = (—4,4)^, £ = 1, a = (1, l)"^, /o = —1 and uq = 0. 
This numerical example is interesting to study as not much is known about blow-up 
problems with non-symmetric spatial operators. The solution behaves as the solution 
to a linear convection-diffusion problem for small t. As time progresses, the nonlinear 
term takes over and the solution begins to exhibit point growth leading to blow-up. 
As in Example 1, we choose to use a small spatial threshold to render the spatial 
contribution to both the error and the estimator negligible. We then reduce the 
temporal threshold and observe how far we can advance towards the blow-up time. 
The results are given in Table [7^ From the data, we conclude that 


l|t^/i||L“(0,r;L“(O)) OC 


Although not much is known about blow-up problems with convection, it is reasonable 
to assume that because the nonlinear term dominates close to the blow-up time that an 
analogous relationship between the magnitude of the exact solution in the L°°{L°°)- 
norm and distance from the blow-up time exists as in Example 1. Assuming that 
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Table 7.2: Example 2 Results 


ttol“*“ 

Time Steps 

Estimator 

Final Time 

\\Uh{T)\\L^^n) 

1 

4 

3.6 

0.78125 

0.886 

0.125 

10 

3.6 

0.97656 

1.322 

(0.125)2 

54 

22.0 

1.31836 

3.269 

(0.125)3 

119 

47.5 

1.41602 

5.107 

(0.125)4 

252 

132.1 

1.48163 

8.059 

(0.125)3 

520 

218.4 

1.51711 

11.819 

(0.125)3 

1064 

664.6 

1.54467 

18.139 

(0.125f 

2158 

1466.1 

1.56224 

27.405 

(0.125)3 

4354 

1421.7 

1.57402 

41.374 

(0.125)9 

8792 

11423.0 

1.58243 

64.450 

(0.125)4° 

17713 

21497.8 

1.58770 

99.190 

(0.125)44 

35580 

21097.1 

1.59092 

145.785 

(0.125)42 

71352 

35862.0 

1.59299 

211.278 


this is indeed the case and following the same reasoning as in Example 1, we again 
conclude that 


A(ttol+,iV) oc 

7.3. Example 3. Let = (—8,8)^, e = 1, a = (0,0)^, /o = 0 and the ‘volcano’ 
type initial condition be given by uq = 10{x^ + y2^g-o.5(£c +y ) ^ blow-up set for 
this example is a circle centred on the origin — this induces layer type phenomena in 
the solution around the blow-up set as the blow-up time is approached making this 
example a good test of the spatial capabilities of the adaptive algorithm. Once more, 
we choose a small spatial threshold so that the spatial contribution to the error and 
the estimator are negligible. We then reduce the temporal threshold and see how far 
we can advance towards the blow-up time. The results are given in Table [73} 


Table 7.3: Example 3 Results 


ttol^" 

Time Steps 

Estimator 

Final Time 

\\UH{T)\\L^(n) 

8 

3 

15 

0.06250 

10.371 

1 

10 

63 

0.09375 

14.194 

0.125 

36 

211 

0.11979 

21.842 

(0.125)2 

86 

533 

0.13412 

31.446 

(0.125)3 

190 

971 

0.14388 

45.122 

(0.125)4 

404 

1358 

0.15072 

64.907 

(0.125)3 

880 

5853 

0.15601 

98.048 

(0.125)3 

1853 

10654 

0.15942 

146.162 

(0.125)^ 

3831 

21301 

0.16176 

219.423 

(0.125)3 

7851 

143989 

0.16336 

332.849 

(0.125)9 

16137 

287420 

0.16442 

505.236 

(0.125)3° 

32846 

331848 

0.16512 

769.652 

(0.125)33 

66442 

626522 

0.16558 

1175.21 
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Fig. 7.2: Example 3: Initial (left) and final (right) meshes. 


Once again, the data implies that 

\\Uh\\L°°{o,T\L°°{n)) oc 
Arguing as in Example 1, we again conclude that 

A(ttol^,iV) (x 


The numerical solution at t = 0 and t = T obtained with the final numerical 
experiment (ttol+ = (0.125)^^) is shown in Figure 7.3 the corresponding meshes are 
displayed in Figure |7.2[ The initial mesh has a relatively homogenous distribution 


of elements which is to be expected since the initial condition is relatively smooth. 
In the final mesh, elements have been added in the vicinity of the blow-up set and 
removed elsewhere, notably near the origin. The distribution of elements in the final 
mesh strongly indicates that the adaptive algorithm is adding and removing elements 
in an efficient manner. 


8. Conclusions. We proposed a framework for space-time adaptivity based on 
rigorous a posteriori bounds for an IMEX dG discretization of a semilinear blow-up 
problem. The error estimator was applied to a number of test problems and appears 
to converge towards the blow-up time in all cases. In Section it was observed 
that the a posteriori error estimator for the related ODE problem with polynomial 
nonlinearity approaches the blow-up time with a rate of at least one for a basic 
Euler method. The numerical examples show that, for the PDE blow-up problem, 
the proposed error estimator appears to be advancing towards the blow-up time at a 
rate approximately half of that observed for the corresponding ODE error estimator. 
A possible reason for this behaviour lies in the proof of the a posteriori bound via 
an energy argument. Nevertheless, it is this very energy argument which delivers 
a practical conditional a posteriori bound in the sense that condition (5.19) can be 
satisfied within a practically relevant (in terms of computational cost) mesh-parameter 
regime. It would be interesting to investigate the derivation of conditional a posteriori 
bounds for fully-discrete schemes for blow-up problems via semigroup techniques, in 
the spirit of [33], although this currently remains a challenging task. 
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Fig. 7.3: Example 3: Initial (above) and final (below) solution profiles. 
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